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ABSTRACT 

A recently developed constitutive model is implemented into LS-DYNA as a user defined 
material model (UMAT) to characterize the nonlinear strain rate dependent behavior of polymers. 
By utilizing this model within a micromechanics technique based on a laminate analogy, an 
algorithm to analyze the strain rate dependent, nonlinear deformation of a fiber reinforced 
polymer matrix composite is then developed as a UMAT to simulate the response of these 
composites under high strain rate impact. The models are designed for shell elements in order to 
ensure computational efficiency. Experimental and numerical stress-strain curves are compared 
for two representative polymers and a representative polymer matrix composite, with the 
analytical model predicting the experimental response reasonably well. 


INTRODUCTION 

In recent years, there is an increasing trend in the use of fiber-reinforced polymer composites to 
primary aircraft structures, such as the radome and the fan containment system of jet engines. 
For all of the above applications, accounting for and correctly computing the strain rate 
dependent mechanical properties of polymer matrix composites are a major concern. 
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To assist in designing composite structures subject to high strain rate impact loads, finite element 
simulation is a powerful tool and can even in some circumstances partially substitute for 
expensive experimental tests. Among the various commercial finite element codes, LS-DYNA [1] 
excels in large deformation transient dynamic problems and impact simulations due to the 
explicit time integration algorithms within the code. 

LS-DYNA has a large library of material options which have been widely used in the automobile 
and aerospace industries. For example, Schweizerhof, et. al. [2] used LS-DYNA to conduct a 
crashworthiness analysis of a composite structure. Bilkhu, et. al. [3] studied the transverse 
impact response of composite plates using LS-DYNA as well. The material models for 
composite materials currently included within LS-DYNA are specifically designed for shell 
elements in order to improve the computational efficiency when conducting a structural analysis 
of large composite structures. However, the material models for composite materials currently 
included within LS-DYNA do not consider the effect of strain rate on the material response, 
which is an important feature in polymer matrix composites, especially under high strain rate 
impact. Recently, Yen [4] implemented a new composite failure model into LS-DYNA. This 
model uses an empirical equation to account for the effect of strain rate on the layer strength 
values within the composite failure criteria, and characterizes the progressive failure behavior. 
However, the effects of strain rate on the deformation response are not accounted for within this 
model, and the model is only applicable to solid elements. 

In this study, a recently developed material model [5] is implemented into LS-DYNA as a user 
defined material (UMAT) designed for shell elements. In this material model, state variable 
constitutive equations originally developed for metals are modified in order to compute the 
nonlinear, strain rate dependent deformation of polymers, including the effects of hydrostatic 
stresses. These equations are then utilized within a strength of materials based micromechanics 
approach using a laminate analogy in order to predict the strain rate dependent deformation of 
polymer matrix composites. Experimentally obtained stress-strain curves of two representative 
polymers and a representative polymer matrix composite under both moderate and high strain 
rate conditions are used to verify the UMAT with a single element simulation. A set of failure 
criteria and a damage progression model are also proposed in this paper. 


POLYMER CONSTITUTIVE EQUATIONS 

To analyze the nonlinear, strain rate dependent deformation of the polymer matrix constituent, 
the Bodner-Partom viscoplastic state variable model [6], which was originally developed to 
analyze the viscoplastic deformation of metals above one-half of the melting temperature, has 
been modified [5,7], In state variable models, a single unified strain variable is defined to 
represent all inelastic strains [8]. Furthermore, in the state variable approach there is no defined 
yield stress. Inelastic strains are assumed to be present at all values of stress, the inelastic strains 
are just assumed to be very small in the “elastic” range of deformation. State variables, which 
evolve with stress and inelastic strain, are defined to represent the average effects of the 
deformation mechanisms. 
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In the modified Bodner model used in this study, the components of the inelastic strain rate 
tensor, el , are defined as a function of the deviatoric stress components S y , the second invariant 

of the deviatoric stress tensor J 2 and an isotropic state variable Z , which represents the 
resistance to molecular flow. The components of the inelastic strain rate are defined as follows 


e 1 = 2 D 0 exp 


/ \ 2n 

' Z ' 




^k +aS ’ 


(i) 


where D 0 and n are both material constants, with D 0 representing the maximum inelastic strain 
rate and n controlling the rate dependence of the material. The effective stress, <j e , is defined as 

(j e =fiT 2 +Sa(J kk ( 2 ) 

where a is a state variable controlling the level of the hydrostatic stress effects and <7 kk is the 
summation of the normal stress components which equals three times the mean stress. 

The rate of evolution of the internal stress state variable Z and the hydrostatic stress effect state 
variable a are defined by the equations 


Z = c/(Z t -Z)e‘ e 

( 3 ) 

a = q (a x - a)e e 

( 4 ) 


where q is a material constant representing the “hardening” rate, and Z, and a x are material 
constants representing the maximum value of Z and a , respectively. The initial values of 
Z and a are defined by the material constants Z 0 and a 0 . The term e[ in Equations 3 and 4 
represents the effective deviatoric inelastic strain rate, which is defined as 


■ i 2 

=Jj e 9 e & 
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where el are the components of the inelastic strain rate tensor and £' m is the mean inelastic 

strain rate. In many state variable constitutive models developed to analyze the behavior of 
metals [8], the total inelastic strain and strain rate are used in the evolution laws and are assumed 
to be equal to their deviatoric values. As discussed by Li and Pan [9], since hydrostatic stresses 
contribute to the inelastic strains in polymers, indicating volumetric effects are present, the mean 
inelastic strain rate cannot be assumed to be zero, as is the case in the inelastic analysis of metals. 
Further information on the constitutive model, along with the procedures required to obtain the 
material constants, can be found in Goldberg, et al [5, 7], 
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COMPOSITE MICROMECHANICAL MODELING 


To compute the effective strain rate dependent, nonlinear deformation response of polymer 
matrix composites based on the response of the individual constituents, a micromechanical 
model has been developed [10]. In this model, the unit cell is defined to consist of a single fiber 
and its surrounding matrix. Due to symmetry, only one-quarter of the unit cell is analyzed. The 
composites are assumed to have a periodic, square fiber packing with a perfect interfacial bond. 
The fibers are assumed to be transversely isotropic and linear elastic with a circular cross- 
section. The matrix is assumed to be isotropic, with a rate dependent, nonlinear deformation 
response computed using the equations described in the previous section of this report. 

The unit cell is divided up into an arbitrary number of rectangular, horizontal slices, as is shown 
in Figure 1. Each slice is assumed to be in a state of plane stress and classical laminate theory is 
assumed to be applicable at this scale. The top and bottom slices in the unit cell are composed of 
pure matrix. The remaining slices are composed of two subslices; one composed of fiber material 
and the other composed of matrix material. For the slices containing both fiber and matrix, the 
out-of-plane stresses can be nonzero in individual subslices, but the volume average of the out- 
of-plane stresses must be equal to zero. By using this approach, the behavior of each slice is 
decoupled, and the response of each slice can be determined independently, which significantly 
reduces the level of complexity in the analysis. Laminate theory is used to obtain the effective 
response of the unit cell. 

The thickness, fiber volume ratio and thickness ratio (the ratio of the slice thickness to the total 
unit cell thickness) for each slice can be determined using the composite fiber volume ratio and 
geometric principles. The overall diameter of the fiber is related to the fiber volume fraction of 
the overall composite, and this term can be used along with the standard geometric definition of 
the radius of a circle to compute the horizontal coordinate of any point on the outer surface of the 
fiber in terms of the fiber volume fraction and the vertical coordinate. Using this information and 
standard calculus concepts, the area of the portion of the fiber contained within each slice of the 
one-quarter of the unit cell which is analyzed can be computed. The fiber diameter and fiber area 
within each slice can be used to compute the fiber volume fraction and thickness ratio of each 
slice. 

The effective properties, effective stresses and effective inelastic strains of each slice are 
computed independently. Micromechanics equations are developed using uniform stress and 
uniform strain assumptions for those slices composed of both fiber and matrix material. The 
stresses and inelastic strains in the slices composed of pure matrix are computed using the matrix 
elastic properties and the inelastic constitutive equations. The standard transversely isotropic 
compliance matrix (or isotropic in the case of the matrix) is used to relate the local strains to the 
local stresses in the fiber and matrix. Each slice is assumed to be in a state of plane stress on the 
global level, but out-of-plane normal stresses can exist in each subslice. The responses of each 
slice are combined using laminate theory to obtain the effective response of the corresponding 
lamina. Laminate theory is applied once again to obtain the effective properties and force 
resultants due to inelastic strains for a multilayered composite laminate. 
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DAMAGE MODEL 


Fiber reinforced composites have a variety of damage mechanisms. At the micro-level, these 
mechanisms can include fiber fracture, fiber buckling, fiber splitting, fiber pullout, fiber/matrix 
debonding and matrix cracking. At the laminate level the damage mechanisms can include 
transverse cracks in planes parallel and perpendicular to the fibers and delamination between 
layers of the laminate [11]. 

Among the various damage modes, Langlie and Cheng [12] stated that fiber breakage, through- 
thickness compressive and shear punching failure and delamination are the major damage 
mechanisms in high strain rate impact. In a high strain rate impact environment, one can argue 
that a strain controlled failure model is more realistic than a stress controlled failure model. 
Therefore, the failure criteria for this study are defined as follows: 


Fiber breakage: 


y: =— si («) 

n 

Compressive punching failure: 

f 2 =~> 1 (?) 

£ f 33 

Shear punching failure: 

/3= _f^ + _fu_> 1 (8) 

^*/23 £ f n 

where e fu is the fiber direction composite tensile failure strain, and £ f3} is the through 

thickness compressive composite failure strain. e f23 and f /13 are the through thickness shear 

composite failure strains. The failure mode of delamination is not included in this material 
model. 

Once failure has occurred, two strategies can be applied to describe the post-failure behavior of 
the material. One option is to assume that the material fails instantaneously, while the other 
option is to apply a post-failure model [13], For the instantaneous failure model, when fiber 
breakage occurs in the 1 direction, all stresses related to that direction (cr. , i = 1 , /= 1 3 ) are set 

to zero instantaneously. When either compressive or shear punching failure occurs, all stresses 
are instantaneously set to zero. 

A post-failure model based on continuum damage mechanics was reported by Matzenmiller, et. 
al. [14], Recent studies conducted using this model show significant improvement in the 


NASA/TM— 2003-212583 


5 



prediction of damage [13], [4], This damage model characterizes the growth of damage by a 
decrease in the stiffness of the material. When fiber breakage occurs in the 1 direction, stiffness 
values in that direction and in all shear directions are reduced. When either compressive or shear 
punching failure occurs, all stiffness values are reduced. It can be described by the following 
equations: 


S = 


(l-^i X 

_Yn_ 

E i 
_Yl3_ 

E \ 

0 


0 
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e 2 
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(i-^K 

E 2 
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( 1-^6 tel 


( 9 ) 


where .S' is the compliance matrix, and the stiffness matrix C is obtained by inverting S, that is 

[Q =[S] _1 

07 1 is the damage variable associated with the z' th failure mode to relate the onset and growth of 
damage to stiffness losses in the material, which can be obtained from an individual failure mode 

j as 


07 , . = 1 - e‘ 




( 10 ) 


where m is the damage exponent, and f } (/=l-3) is the damage parameter as described in 
equations (6)-(8). At the onset of failure, the value of 07 i is zero. As damage progresses, 07 : 

increases and therefore eventually decreases the stiffness of the material until it reaches a final 
value of zero. 


LS-DYNA IMPLEMENTATION AND MODEL VERIFICATION 

LS-DYNA has a user defined material option where a user can implement his or her own 
material model into the code through the use of a UMAT subroutine. The structure of the 
subroutine is shown in the Figure 2. The LS-DYNA code calculates the strain increments for a 
time step and passes them to the UMAT subroutine at the beginning of each time step. The 
material constants, such as moduli and Poisson’s ratios, are read from the LS-DYNA input file 
by the subroutine. The current value of the state variables such as Z and a are saved as history 
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variables which can be read in by the subroutine. By using the history variables, material 
constants and strain increments, the subroutine is able to calculate the stresses at the end of the 
time step by using an incremental form of the polymer constitutive equations and the composite 
micromechanics equations. The subroutine then updates and saves the history variables for use in 
the next time step, and outputs the calculated stresses together with the through thickness strain 
increment. For shell elements under plane stress conditions, the through thickness strain 
increment for the time step must be computed in the UMAT subroutine. The subroutine also 
checks the damage and modifies the material stiffness if damage occurs. The procedure 
described above is carried out on each integration point of each element independently. 

As a first step in implementing the analytical model into LS-DYNA, the polymer constitutive 
model was implemented into the code as a UMAT. To verify the analytical model and its 
implementation within LS-DYNA, two representative toughened epoxies, PR520 and 977-2, 
were analyzed using LS-DYNA, and the computed stress-strain curves were compared to 
experimental results obtained by Goldberg, et. al. [5, 7]. In the experimental tests, longitudinal 
tensile tests and pure shear tests were conducted at room temperature on the materials at strain 
rates of about 5xl0 -5 /sec, 1 /sec and 400 /sec. The low and moderate strain rate tests were 
conducted using an Instron hydraulic testing machine. The high strain rate tests were conducted 
using a split Hopkinson bar. Engineering stress and engineering strain were measured until 
failure. 

In the finite element analyses, a four-node single Belytschko-Tsay shell element was used for the 
simulation of the pure shear and tension tests. The load and boundary conditions applied to the 
model are shown in Figure 3. The applied strain rate is controlled by the displacement-time curve 
defined in the LS-DYNA input file. Since LS-DYNA is an explicit code, designed for solving 
dynamic problems, attempts to simulate the extremely low stain rate tests were found to have 
numerical stability problems, and the analyses took an extremely long time to reach the failure 
strain, even by using only a single element. In addition, in a realistic impact problem, the strain 
rates are well above the quasi-static level. Therefore, the low strain rate case of 5xl0 -5 /sec was 
not considered here. Table 1 lists the material constants used in the model, which were 
determined using the procedures described in Goldberg, et. al. [5, 7]. 


TABLE 1.— MATERIAL CONSTANTS USED IN POLYMER MATRIX MODELING 



Strain 

Rate 

(1/sec) 

Modulus 

(GPa) 

Poisson’s 

Ratio 

D 0 

(1/sec) 

n 

Zo 

(MPa) 

Zi 

(MPa) 

q 

a 0 

OCi 

PR520 

1.76 

3.54 

0.38 

lxlO 6 

0.93 

396.09 

753.82 

279.26 

0.568 

0.126 

420 

7.18 

977-2 

1.9 

3.52 

0.40 

lxlO 6 

0.85 

259.50 

1131.4 

150.50 

0.129 

0.152 

500 

6.33 


Experimental [5,7] and computed shear stress-shear strain curves for PR520 are shown in Figure 
4 for the moderate and high strain rate cases, while tensile stress-strain curves are shown in 
Figure 5. Experimental [5] and computed shear stress-shear strain curves for 977-2 are shown in 
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Figure 6, and tensile stress-strain curves are shown in Figure 7. Both materials exhibit a strain 
rate dependent, nonlinear deformation response under both shear and tensile loading. For the 
experimental data, at high strain rates, the sharp increase in stress at the beginning of the loading 
with negligible increase in strain observed for PR520 under shear and tensile loading is most 
likely not representative of the actual material behavior, but rather an artifact of the experimental 
tests. Also note that the oscillations observed in the high strain rate tensile tests for 977-2 were a 
result of the specimen geometry that was utilized for these tests. 

From Figures 4-7, it can be seen that the finite element results under shear loading correlate 
reasonably well with the experimental data for both materials. The nonlinearity and rate 
dependence observed in the deformation of the two polymers are captured qualitatively, and the 
quantitative match between the experimental and computed results is reasonably good. For the 
tensile loading condition, while the strain rate dependence and nonlinearity of the stress-strain 
results are correctly captured qualitatively, for the PR520 material at the moderate strain rate the 
stresses in the nonlinear portion of the stress-strain curve are under predicted for reasons that 
require further investigation. One possible cause for the discrepancy may be that the beginning 
and final values of the hydrostatic stress state variable a (ao and oq) actually vary with strain rate, 
while in this study these values were assumed to be constant and were determined using low 
strain rate data not included here. For the PR520 material under high strain rate tensile loading, 
the slope of the nonlinear portion of the stress-strain curve matches the slope of the experimental 
curve, and if the experimental curve was shifted appropriately the comparison between the 
experimental and computed results may actually be quite good. For the 977-2 material under 
tensile loading, at the moderate strain rate in the nonlinear portion of the stress-strain curve the 
stresses are still slightly under predicted, but the comparison between the experimental and 
computed results are much improved as compared to the PR520 predictions. For 977-2 under 
high strain rate tensile loading, the predicted stress-strain curve bisects the experimental curve, 
indicating that the predictions may adequately represent the actual material behavior. 

The next step in the analysis process involved implementing the composite micromechanical 
model into LS-DYNA through the use of a UMAT user defined subroutine. A series of analysis 
studies were then carried out to verify the finite element implementation of this model. 

The composite material used for the verification studies consists of carbon IM7 fibers in the 
977-2 toughened epoxy matrix discussed earlier. Each composite was defined to have four plies 
with a fiber volume of 0.60. The properties for 977-2 matrix are the same as defined in Table 1. 
The IM7 fiber properties are as follows: En=276 GPa, E 2 2=13.8 GPa, Vi 2 = 0 . 25 , V 23 = 0 . 25 , 
Gi2=20.0 GPa. The loading and boundary conditions used for the composite verification studies 
are identical to those used for the verification studies for the polymer constitutive equations. 

Experimental [5] and predicted longitudinal tensile stress-strain curves for two laminate 
configurations ([45°] and [90°]) of the IM7/977-2 material at strain rates of approximately 1 /sec 
and 400 /sec are shown in Figure 8 and Figure 9. In Figure 8, for the [45°] laminate at the high 
strain rate the LS-DYNA stress-strain curve exhibited a wave shape, which may due to the shock 
wave from the boundary induced by the high strain rate load on the single element. This 
discrepancy will be investigated further. In Figure 9, the experimental curve shows a higher 
stress at the initial stage of loading than the numerical results. However, the high stress levels 
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seen in the experimental curves may be an artifact of the way the stress waves progress through 
the specimen at the start of the experiment, and the computations may actually be fairly accurate. 
At higher strain levels, the comparison is much improved. The comparison between the 
experimental and numerical results is excellent for the moderate strain rate in both of the 
laminate configurations. 

A preliminary attempt to exercise the damage models is shown in Figure 10 for the progressive 
damage model. The bulk polymer PR520 was analyzed under tensile loading at the moderate 
strain rate of 1.4 /sec. It can be seen that by adjusting the damage exponent m, damage 
progression can be controlled. Large values of m represent brittle fracture, whereas lowering the 
value of m increases the post- failure strength. 


CONCLUSIONS 

State variable constitutive equations based on the Bodner viscoplasticity model have been 
modified to analyze the nonlinear, strain rate dependent deformation of polymers. The effects of 
hydrostatic stresses on the inelastic deformation have been accounted for. The constitutive 
equations have been implemented within the LS-DYNA transient dynamic finite element code 
through the use of user defined subroutines (UMATs). The tensile and shear deformation of two 
representative polymers have been accurately simulated using the polymer UMAT with the 
constitutive model. 

The polymer constitutive equations were then implemented into a composite micromechanical 
model where the unit cell has been divided into several independently analyzed slices. The 
composite model has also been implemented into LS-DYNA through the use of a UMAT. The 
tensile deformation of a representative polymer matrix composite was successfully predicted for 
two fiber orientations and two strain rates, indicating that the composite UMAT is able to capture 
the important features of the deformation response. There were some discrepancies between the 
experimental and predicted responses when high strain rate tests were simulated. The causes of 
these discrepancies will be investigated further. 

A set of damage criteria and a damage progression model was also proposed in this paper. The 
damage progression model was implemented into the polymer UMAT as an initial attempt to 
develop a damage model for the entire composite. Once the damage model is developed to the 
point where it can be applied to the constituents of the composite, the damage model has great 
potential to predict and simulate the composite damage progression. Therefore, the damage 
progression model will be further developed and implemented into the composite UMAT. 

The current composite UMAT is based on a plane stress assumption. Future work will 
incorporate the out-of-plane transverse shear stresses in order to more accurately predict the 
material deformation under impact conditions. Since the composite UMAT is based on 
micromechanics, the model will be modified to allow for the analysis of woven and braided 
composites. 
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Figure 1. — Unit cell and slices in the micromechanical model 
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Figure 2. — Flowchart of user defined subroutine used to implement material 
model within LS-DYNA finite element code. 
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Figure 3. — Boundary and loading conditions for single element 
(a) shear and (b) tension finite element analyses 



Strain 


Figure 4. — Experimental and LS-DYNA simulated shear stress-shear strain curves 
for PR520 resin at strain rates of 1.76 /sec and 420 /sec 
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Figure 5. — Experimental and LS-DYNA simulated tensile stress-strain curves 
for PR520 resin at strain rates of 1.4 /sec and 510 /sec 
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Ls-dyna: 1.91/s 
Ls-dyna: 518/s 
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Figure 6. — Experimental and LS-DYNA simulated shear stress-shear strain curves 
for 977-2 resin at strain rates of 1.91 /sec and 518 /sec 
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Figure 7. — Experimental and LS-DYNA simulated tensile stress-strain curves 
for 977-2 resin at strain rates of 1.31 /sec and 365 /sec 



+ Experiment: 1.2/s 
Ls-dyna: 1.2/s 
a Experiment: 405/s 
Ls-dyna: 405/s 


Figure 8. — Experimental and LS-DYNA simulated tensile stress-strain curves 
for IM7/977-2 [45°] composite at strain rates of 1.2 /sec and 405 /sec. 
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Strain 


Figure 9. — Experimental and LS-DYNA simulated tensile stress-strain curves 
for IM7/977-2 [90°] composite at strain rates of 1.09 /sec and 405 /sec. 



Strain 

Figure 10. — Examination of the effects of damage exponent m on the damage progression 
for PR520 resin under tensile loading at a strain rate of 1 .4 /sec. 
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